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Abstract 

A quantum equation-of-state model is presented and applied to the 
calculation of high-pressure shock Hugoniot curves beyond the asymp- 
totic fourfold density, close to the maximum compression where quan- 
tum effects play a role. An analytical estimate for the maximum at- 
tainable compression is proposed. It gives a good agreement with the 
equation-of-state model. 

1 Introduction 

Maximum pressures reached nowadays by the shock-wave technique in labo- 
ratory laser experiments are as high as hundreds of megabars and even more 
for many materials. As the strength of the shock is varied for a fixed ini- 
tial state, the pressure-density final states of the material behind the shock 
belong to a curve named shock adiabat or Hugoniot curve. The maximum 
compression attainable by a single shock is finite and occurs for a finite 
pressure. This phenomenon is due to the draining of energy in internal de- 
grees of freedom, via ionization and excitation. In this region, the electrons 
from the ionic cores are being ionized and the shock density increases be- 
yond the fourfold density 4po (po being inital density), which corresponds 
to an asymptotic situation where ionization is completed and the plasma 
approaches an ideal gas of nuclei and electrons. The Hugoniot curve, and 
therefore the maximum compression, depend on the equation of state (EOS) 
of the matter, which, in principle, can be determined from theory. In prac- 
tice, the barriers to ah initio calculations are formidable owing to the com- 
putational difficulty of solving the many-body problem. Consequently, it has 
proven necessary to introduce simplifying approximations into the govern- 
ing equations. Many EOS models rely on the Thomas-Fermi (TF) model of 
dense matter [Q. This approach contains certain essential features in order 
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to characterize the material properties but is valid over a limited range of 
conditions and provides only a semi-classical description of electrons. How- 
ever, at intermediate shock pressures, when the material becomes partially 
ionized, the EOS depends on the precise quantum-mechanical state of the 
matter, i.e. on the electronic shell structure. 

In Sec. El a quantum self-consistent-field (QSCF) EOS model is pre- 
sented. The potential and the bound-electron wavefunctions are determined 
solving Schrodinger equation by a self-consistent procedure in the density 
functional theory using a finite-temperature exchange-correlation potential. 
Hugoniot curves are calculated and focus is put on the region where quan- 
tum shell effets are visible, in the vicinity of the maximum compression. In 
Sec. El an analytical expression for the maximum compression attainable by 
a single shock is proposed. It can be applied to any material from any initial 
state except those with gaseous densities. Maximum compression evaluated 
by this formula is compared to the one obtained from our QSCF model 
presented in Sec. El 

2 Quantum self-consistent-field equation-of-state 
model 

Hugoniot curves are obtained by resolution of equation 

(P + P )(l/p -l/p) = 2(E-E ), (1) 

which requires the knowledge of pressure P and internal energy E versus 
density p and temperature T over a wide region of the phase diagram. Sub- 
script characterizes the initial state. In the present work, focus is put on 
the standard Hugoniot (Pq=0, Po equal to solid density and To=300 K). The 
adiabatic approximation is used to separate the thermodynamic functions 
into electronic and ionic components. The total pressure of the plasma can 
be written 

P(p,T) = P e (p,T) + P i (p,T), (2) 

where P e (p, T) and Pi(p, T) are respectively the electronic and ionic con- 
tributions to the pressure. The electronic contribution to the EOS stems 
from the excitations of the electrons due to temperature and compression. 
Atoms in a plasma can be idealized by an average atom confined in a Wigner- 
Seitz (WS) sphere, which radius r ws is given by solving 

gTTw^-^Oo = !> ( 3 ) 

where ao=52. 9177208319 10 -8 cm is the Bohr radius. Na the Avogadro 
number, p the matter density in g/cm 3 , and A the atomic weight in g. 
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Throughout the paper, all other quantities are in atomic units (a.u.) defined 
by m = h = e = 1, m and e being respectively electron mass and charge. 
Inside the sphere, the electron density has the following form [2]: 



n(r) = fiitb ~ M) ^(f) 2 + ^ {kB p 3/2 J 1/2 ( - V(r),p - V(rj) (4) 



b 

where 



1 + exp|x/ kbT\ 
is the usual Fermi-Dirac population and 

J n /2{a,x)= — -dy (6) 

' J a 1 + exp(y - x) 

is the incomplete Fermi function of order n/2. The first term in (j3J) is 
the contribution of bound electrons to the charge density, while the second 
term is the free-electron contribution, written in its semi-classical TF form. 
The wavefunction ifi b associated to energy e& of a bound orbital is solution 
of the Schrodinger equation 

- -Aij} b + V{r)ip b = e b i; b , (7) 



V(r) = ksTV(r) being the self-consistent potential: 

V(r) = --+ -^L*V + V xc (r), (8) 



r 



where V xc is the exchange-correlation contribution, evaluated in the local 
density approximation [3j. Last, the chemical potential \i is obtained from 
the neutrality of the ion sphere: 

IS 

n(r)Airr 2 dr = Z, (9) 

and p = n/(k B T). Eqs. g]), (® and © must be solved self- 
consistently. The electronic pressure P e [U[5] consists of three contributions, 
P e = P b + Pf + P xc , where the bound-electron pressure P b is evaluated at 
the boundary of the WS sphere using the stress-tensor formula 



p \- fi(e b -n) \(dy b \ 2 ,f l + l + l 2 \ 2f 1 

yb being the radial part of the wavefunction ip b multiplied by r. The 
free-electron pressure Pf reads 
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P f = ^(^ T ) 5/2j 3/ 2 ( " V(r ws ),Jl- V(r ws )) 

and P xc is the exchange-correlation pressure evaluated in the local den- 
sity approximation [3 J . The value of bound-electron pressure P& depends on 
the boundary conditions. This comes from the fact that the energy of an 
orbital depends on the value of the corresponding wavefunction and on its 
derivative at the boundary of the WS sphere. Tabled] represents the energies 
of the orbitals of iron (Fe, Z=26) obtained if the wavefunction cancels at the 
boundary (BC1), if the logarithmic derivative of the wavefunction cancels 
at the boundary (BC2) and if the wavefunction behaves like a decreasing 
exponential at the boundary (BC3). In this example, the mass density is 
equal to 34.83 g/cm 3 and the temperature to 115.30 eV. Such conditions 
belong to the Hugoniot of Fe (see Fig. [2]) calculated from the present QSCF 
model. The values of orbital energies displayed in table [1] show that the 
orbitals calculated with the three different boundary conditions will not be 
ionized for the same value of density. We chose to perform single-state cal- 
culations with boundary condition (BC3), since it allows the matching of 
wavefunctions with Bessel functions outside the WS sphere [I], where the 
potential is zero. Moreover, condition (BC3) is consistent with the fact that 
the bound-electron density vanishes at infinity. The internal energy is 

Ee = S- nk e k -\ n(r) ^X-dW 
Z ^ Jo Jo \r — r'\ 

rr wa 

+E XC - / n(r)V xc (n(r))d 3 r, (11) 
Jo 

where E xc is the exchange-correlation internal energy and n k the popula- 
tion of state k (either bound or free). The first term in (jlip can be expressed 
by 

^n k e k = ^2 fi(e b , /JL)e b + E fik + E f>p . (12) 

k b 

Ef :P is the free-electron potential energy 

Ef,P = ^ ( y 3/2 [J S Ji/2 ( " V(r),H- V(r))v(r)d 3 r (13) 
and Ef tk the free-electron kinetic energy 

= V2(k B p^ r» j(_ - (r)) _ _ _ (r) x 3 r (M) 

^ Jo 7 
The ionic contribution to the EOS can be estimated in the ideal-gas 
approximation including non-ideality corrections to the thermal motion of 
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ions. We used an approximation based on the calculation of the EOS of 
a one-component plasma (OCP) by the Monte Carlo method [6]. The ion 
contribution can be obtained using a formula based on the Virial theorem: 

Pi(j>, T) = pk B T + ^AEiip, T), (15) 

where 

AEkQ,, T) = k B T [F 3 / 2 £ - + a - r)t/2 - ai r] (16) 

and Ei(p,T) = 3k B T/2+AEi(p,T) withoi = -0.895929, a 2 = 0.11340656, 
a 3 = -0.90872827, a 4 = 0.11614773, h = 4.666486, b 2 = 13.675411, 
63 = 1.8905603 and 64 = 1.0277554. T is the plasma coupling parameter 
(ratio of ionic Coulomb interaction and thermal kinetic energy). In the 
present study, T > 1. 

The contributions P e (p, T) and E e (p, T) calculated from the QSCF model 
are not valid for T — > 0. Therefore, P e (p, 0) and E e (p, 0) must be substracted 
from total pressure and internal energy. They are replaced by P c (p) and 
E c (p), which constitute the cold curve. The total pressure reads now 

P(p,T) = P e (p,T) + P(p,T) - P e (p,0) + P c (p). (17) 

In our model, the cold curve is obtained either from Augmented Plane 
Waves (APW) [7] simulations, or using the Vinet [8] universal EOS. In many 
situations, Vinet EOS gives realistic results, and its accuracy is sufficient for 
most applications involving high pressure and high temperature situations. 

Fig. [T] represents the principal Hugoniot curve for Al calculated from a 
model relying on Thomas-Fermi approximation and from the present QSCF 
model. Fig. [2] illustrates the same calculations in the case of Fe. One can 
check that the maximum compression is beyond the ideal-gas asymptote 
p = 4po- The oscillations are a consequence of the competition between the 
release of energy stocked as internal energy within the shells and the free- 
electron pressure. When ionization begins, the energy of the shock is used 
mainly to depopulate the relevant shells and the material is very compres- 
sive. However, the pressure of free electrons in increasing number dominates 
again and the material becomes more difficult to compress. The two shoul- 
ders in the Hugoniot of Al are due to the successive ionization of L and K 
shells respectively. Similarly, the three shoulders in the Hugoniot of Fe are 
due to the successive ionization of M, L and K shells respectively. Table [2] 
illustrates the impact of the choice of the boundary conditions of the wave- 
functions and of the calculation of the ionic part (perfect ideal gas (IG) with 
or without OCP corrections (|15til6p ) on the maximum compression rate. As 
for the eigen-energies (see table [[]), (BC1) gives the highest value, (BC2) 
the lowest and (BC3) the intermediate one. The differences concerning the 
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maximum compression rate appear to be less than 1 %. OCP corrections 
systematically increase the maximum compression rate. 



3 Simple evaluation of maximum compression 

The total energy can be written as the sum of kinetic and potential ener- 
gies E k and E p . Neglecting exchange-correlation contribution for a sake of 
simplicity, the virial theorem enables one to relate pressure, kinetic energy 
and potential energy: 

3- = 2E k + E p . (18) 
P 

Using the Hugoniot relation ([I]), the compression rate rj = p/po for the 
standard Hugoniot (Pq = 0, po solid density and To=300 K) can be written 

3 3 
? ? = 4 + =4 + (19) 

l + 2 E^E^ L + Z AE P 

At high compression, assuming that E k » E ko and considering that all 
the electrons have been ionized and that all electrons have a kinetic energy 
equal to the Fermi energy, it is possible to write: 

^E k = z\{^Zp^)^al (20) 

At high compression, the excess potential energy can be estimated as 
the Coulomb interaction energy of two ionic spheres at close contact 



1 Z 



2 



3A 



1/3 



p- 1/3 , (21) 



AE n = and r„ 

where r ws is the WS radius defined by Eq. ©. Eqs. (f2"U|) and (f2"T|) 
are relevant for a strongly coupled gas of degenerate electrons. Therefore, 
putting Eqs. (|20p and (|2ip in Eq. (|19p . the maximum compression rate rj m 
obeys the following equation 



4 + 573, (22) 

1 + 'y(j> ,Z, A)r]rA 



with 



7 (po,Z,A) =3vr(2iV A ) 1 /3 ao (|^) 1/3 . (23 ) 
Intermediate variable X = r]H 3 obeys the following fourth order equation 
7X 4 + X 3 - 4 7 X - 7 = 0. (24) 
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The solution is 



Vm = 




47 3 \/M7) J 



32 7 3 - 1 -1 3 



(25) 



with e 



1 if 7 < 0.314980262473 and e = +1 else, 



1 2 10 /3 A x /3( 7 ) 



(26) 



4 7 2 AV3( 7 ) ^ 2V3 7 



and 



A( 7 ) = -7 + I67 3 + ^49 + 1824 7 3 + 256 7 6 . 



(27) 



Neglecting cohesive and dissociation energies, and using fits for the total 
ionization energies, J.D. Johnson [9] has proposed an analytical formula for 
the maximum compression rate. Fig. displays the maximum compres- 
sion rate obtained from our analytical formula, from Ref. [9] and calculated 
with our QSCF model for Be, B, C, Na, Mg, Al, Fe and Cu. Exact val- 
ues for those particular elements are indicated in table El It first confirms 
the fact that the maximum compression is always smaller than 7 [9], and 
strongly dependent on the density pq. The most important point is that 
the maximum compression obtained with formula (|25p is very close to the 
one obtained from our QSCF model. This means that Eq. (|25p . which does 
not rely on quantum mechanics, gives a good agreement with a quantum 
EOS model. At first sight, it seems difficult to decide in an unequivocal way 
which approach, between formula (|25p and Ref. [9], is the most accurate. 
However, our estimate (|25p does not rely on a particular form of the EOS 
and is fully analytical, since it does not require any data concerning the 
ionization energies, as in [§]. On the contrary to the formula proposed in 
[9], the maximum compression predicted by our analytical model is higher 
for Fe than for Al, which is consistent with the results presented in |10j . 
However, in [10] the authors notice that the maximum compression seems 
to increase with Z, even if they confess that they have no explanation for 
that. We believe this is only a global tendency, and Fig. [3] shows that the 
maximum compression, evaluated either from Ref. [9] or from Eq. (|25p . 
does not vary monotonically with atomic number Z. 

4 Conclusion 

A quantum equation-of-state model was presented. It consists in a self- 
consistent calculation of the electronic structure. Bound electrons are treated 
in the framework of quantum mechanics and bound-electron pressure is eval- 
uated using the stress-tensor formula. Free electrons are described in the 
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Thomas-Fermi approximation. Exchange-correlation effects at finite tem- 
perature are taken into account. The ionic part is described in the ideal- 
gas approximation with non-ideality corrections from the one-component 
plasma. It was shown that such a model is well suited for the computation 
of Hugoniot curves and exhibits oscillations due to the ionization of suc- 
cessive shells. A simple analytical estimate was proposed for the maximum 
compression attainable by a single shock, which is crucial for the diagnostic 
of high-pressure laser-induced shock waves. It relies on an expression for the 
increase of kinetic energy assuming degenerate electrons and the increase of 
potential energy is evaluated as Coulomb interation energy between hard 
spheres. Furthermore, the maximum compression obtained from this for- 
mula is in good agreement with the one calculated from our quantum EOS 
model. 
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Orbital 


e^'i) ( a .u.) 


(a.u.) 


(a.u.) 


Is 


-256.72931 


-257.32658 


-256.98981 


2s 


-29.891050 


-30.363675 


-30.094431 


2p 


-25.511960 


-26.002778 


-25.723564 


3s 


-2.3491653 


-2.8962546 


-2.5844217 


3p 


-1.0321504 


-1.6990494 


-1.3273471 



Table 1: Energies of orbitals of Fe at /3=34.83 g/cm 3 and T=115.30 eV for 
three boundary conditions: (BC1): the wavefunction is zero at the bound- 
ary, (BC2): the logarithmic derivative of the wavefunction is zero at the 
boundary and (BC3): the wavefunction behaves like a decreasing exponen- 
tial at the boundary. 



Boundary cond. 


Ionic part 


Max. compression rate 


(BC1) 


IG 


4.908 


(BC2) 


IG 


4.891 


(BC3) 


IG 


4.901 


(BC3) 


IG+OCP corr. 


4.931 



Table 2: Impact of the boundary conditions of the bound-electron wave- 
functions and of the treatment of ions on the maximum compression for 
Al. 



Element 


Z 


r] m [QSCF) 


rj m [Formula (125ft] 


WRef. [9J] 


Beryllium 


4 


4.56 


4.70 


4.31 


Boron 


5 


4.61 


4.73 


4.36 


Carbon 


6 


4.60 


4.76 


4.40 


Sodium 


11 


5.07 


5.23 


5.63 


Magnesium 


12 


5.20 


5.13 


5.27 


Aluminum 


13 


4.90 


5.07 


5.06 


Iron 


26 


5.12 


5.15 


5.02 


Copper 


29 


5.14 


5.18 


5.04 



Table 3: Maximum compression rate for different elements calculated from 
our QSCF equation-of-state model, from our analytical formula (|25p and 
from an estimate published in Ref. [S]. 
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Figure 1: Standard Hugoniot curve for Al (po =2.70 g/cm 3 ) calculated from 
a model relying on Thomas-Fermi approximation and from our QSCF EOS 
model. 




Compression 

Figure 2: Standard Hugoniot curve for Fe (po =7.85 g/cm 3 ) calculated from 
a model relying on Thomas-Fermi approximation and from our QSCF EOS 
model. 
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Atomic number Z 

Figure 3: Maximum compression rate obtained from [9], from our model and 
compared to the maximum compression calculated from our QSCF model 
for Be, B, C, Na, Mg, Al, Fe and Cu. 
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